# Load data
data <- read_csv("Data/meta_analysis/existing_tests.csv", show_col_types = FALSE)

# Fill NAs in study variable
data <- data %>%
    mutate(
        study = study %>% replace_na("NA")
    )

# Count number of unique tests, studies, and papers
## Tests
data %>%
    nrow()

## Studies
data %>%
    select(authors, study) %>%
    unique() %>%
    nrow()

## Papers
data$authors %>%
    unique() %>%
    length()

# Lump treatments
data <- data %>%
    mutate(
        `Type of Correction` = case_when(
            treatment_type == "Policy opposition + Meta-dehumanization" ~ "Misc.",
            treatment_type == "Policy opposition" ~ "Policies",
            grepl("dehuman|affect", treatment_type) ~ "Animosity",
            grepl("undemocratic", treatment_type) ~ "Undemocratic Practices",
            grepl("violence", treatment_type) ~ "Violence",
            treatment_type == "Demographics" ~ "Demographics"
        )
    )

# Lump DVs
data <- data %>%
    mutate(
        `Type of Outcome` = case_when(
            dv_type == "Social distance" ~ "Partisan Animosity",
            dv_type == "Partisan animosity" ~ "Partisan Animosity",
            dv_type == "Support for undemocratic practices" ~ "Support for Undemocratic Practices",
            dv_type == "Support for partisan violence" ~ "Support for Partisan Violence"
        )
    )

# Flip betas
data <- data %>%
    mutate(
        flipped_beta = case_when(
            proper_coef_sign == "positive" ~ beta,
            proper_coef_sign == "negative" ~ beta * -1
        )
    )

# Standardize variables
data$standard_beta <- data$flipped_beta / data$dv_range
data$standard_se <- data$se / data$dv_range

# Generate "significance" column
data <- data %>%
    mutate(
        Significance = (abs(beta) / se) > 1.96,
        Significance = Significance %>%
            factor(
                labels = c(
                    "p > .05",
                    "p < .05"
                )
            ),
        significant_and_correct = (flipped_beta / se) > 1.96
    )

# Calculate facet-wise parameters
facet_wise_data <- data %>%
    group_by(`Type of Outcome`) %>%
    reframe(
        max_beta = max(standard_beta, na.rm = TRUE),
        standard_beta = weighted.mean(
            x = standard_beta,
            y = 1 / standard_se,
            na.rm = TRUE
        ),
        percent_significant_and_correct = mean(significant_and_correct, na.rm = TRUE) * 100
    ) %>%
    mutate(
        y_intercept_left = -1 / 1.96 * -standard_beta,
        y_intercept_right = 1 / 1.96 * -standard_beta
    )

####################
### FUNNEL PLOTS ###
####################

# Set visual parameters
source("Code/ggplot2_theme.R")
prl_palette <- c("#E69F00", "#56B4E9", "#009E73", "#0072B2", "#D55E00", "#CC79A7")

# Plot
data %>%
    ggplot() +
    geom_section(
        data = facet_wise_data,
        aes(
            slope = -1 / 1.96,
            intercept = y_intercept_left * 100,
            above = TRUE
        )
    ) +
    geom_section(
        data = facet_wise_data,
        aes(
            slope = 1 / 1.96,
            intercept = y_intercept_right * 100,
            above = TRUE
        )
    ) +
    geom_point(
        data = data,
        aes(
            x = standard_beta * 100,
            y = standard_se * 100,
            color = `Type of Correction`,
            shape = Significance
        ),
        size = 4
    ) +
    geom_blank(
        data = data,
        aes(
            x = standard_beta * 100 * 1.03,
            y = standard_se * 100 * 1.1
        )
    ) +
    geom_vline(
        data = facet_wise_data,
        aes(xintercept = standard_beta * 100),
        linetype = "dotted"
    ) +
    geom_abline(
        data = facet_wise_data,
        aes(
            intercept = y_intercept_left * 100,
            slope = -1 / 1.96,
            group = `Type of Outcome`
        ),
        linetype = "solid"
    ) +
    geom_abline(
        data = facet_wise_data,
        aes(
            intercept = y_intercept_right * 100,
            slope = 1 / 1.96,
            group = `Type of Outcome`
        ),
        linetype = "solid"
    ) +
    geom_text(
        data = facet_wise_data,
        aes(
            label = paste0(
                "p < .05: ",
                round(
                    percent_significant_and_correct,
                    1
                ),
                "%"
            )
        ),
        x = Inf,
        y = Inf,
        hjust = 1.02,
        vjust = 1.3
    ) +
    scale_shape_manual(
        values = c(10, 16)
    ) +
    scale_color_manual(
        values = prl_palette
    ) +
    # scale_x_continuous(limits = c(-1, NA), expand = c(0, 0)) +
    scale_y_continuous(
        limits = c(NA, 0),
        expand = c(0, 0),
        transform = "reverse",
        labels = number_format(accuracy = 0.1)
    ) +
    facet_wrap(
        ~`Type of Outcome`,
        nrow = 3,
        scales = "free_y"
    ) +
    labs(
        x = "**Effect in Desired Direction**<br>(Beta as % of Outcome's Range)",
        y = "**Study Precision**<br>(SE as % of Outcome's Range)<br>",
    ) +
    guides(
        shape = guide_legend(
            nrow = 2,
            byrow = TRUE,
            override.aes = list(
                size = 4
            ),
            title.position = "top"
        ),
        color = guide_legend(
            nrow = 2,
            byrow = TRUE,
            override.aes = list(
                shape = 15,
                size = 4
            ),
            title.position = "top"
        )
    ) +
    theme_prl() +
    theme(
        axis.title.x = element_markdown(),
        axis.title.y = element_markdown(),
        legend.position = "bottom",
        legend.box = "horizontal",
        legend.justification = "center",
        legend.text = element_text(size = 9.5)
    )

## Save plots
ggsave(
    "Images/funnel_plots.png",
    width = 6, height = 9
)

######################
### CORRESPONDENCE ###
######################

# Check means
data %>%
    group_by(treat_dv_identical) %>%
    reframe(
        standard_beta = mean(standard_beta * 100)
    )

# Run model
correspondence_model <- data %>%
    mutate(standard_beta_100 = standard_beta * 100) %>%
    lm(standard_beta_100 ~ treat_dv_identical, data = .)

correspondence_model %>%
    summary()

# Regression table
modelsummary(
    correspondence_model,
    coef_rename = c(
        "(Intercept)",
        "Correction and Hostile Attitude are Near-Identical"
    ),
    title = "Effect of Corrections on Matching versus Non-Matching Hostile Attitudes",
    stars = TRUE,
    output = "Tables/meta_analysis_correspondence_effect.txt"
)

#############
### TABLE ###
#############

data %>%
    rename(
        Authors = authors,
        Study = study
    ) %>%
    arrange(Authors, Study, `Type of Outcome`, `Type of Correction`) %>%
    select(
        Authors, Study, `Type of Outcome`, `Type of Correction`, standard_beta,
        standard_se
    ) %>%
    mutate(
        standard_beta = standard_beta * 100,
        standard_se = standard_se * 100
    ) %>%
    rename(
        "Beta (in Desired Direction)" = standard_beta,
        "SE" = standard_se
    ) %>%
    datasummary_df(
        title = "Tests in Meta-Analysis",
        notes = "Note: Betas and standard errors were standardized by dividing
        these quantities by the outcome variable's range and multiplying by 100.",
        output = "Tables/tests_in_meta_analysis.txt"
    )
